Coursework 1: Image filtering¶

In this coursework you will practice techniques for image filtering. The coursework includes coding questions and written questions. Please read both the text and the code in this notebook to get an idea what you are expected to implement.

What to do?¶

  • Complete and run the code using jupyter-lab or jupyter-notebook to get the results.

  • Export (File | Save and Export Notebook As...) the notebook as a PDF file, which contains your code, results and answers, and upload the PDF file onto Scientia.

  • Instead of clicking the Export button, you can also run the following command instead: jupyter nbconvert coursework_01_solution.ipynb --to pdf

  • If Jupyter complains about some problems in exporting, it is likely that pandoc (https://pandoc.org/installing.html) or latex is not installed, or their paths have not been included. You can install the relevant libraries and retry. Alternatively, use the Print function of your browser to export the PDF file.

  • If Jupyter-lab does not work for you at the end (we hope not), you can use Google Colab to write the code and export the PDF file.

Dependencies:¶

You need to install Jupyter-Lab (https://jupyterlab.readthedocs.io/en/stable/getting_started/installation.html) and other libraries used in this coursework, such as by running the command: pip3 install [package_name]

In [1]:
# Import libaries (provided)
import imageio.v3 as imageio
import numpy as np
import matplotlib.pyplot as plt
import noise
import scipy
import scipy.signal
import math
import time

1. Moving average filter (20 points).¶

Read the provided input image, add noise to the image and design a moving average filter for denoising.

You are expected to design the kernel of the filter and then perform 2D image filtering using the function scipy.signal.convolve2d().

In [2]:
# Read the image (provided)
image = imageio.imread('campus_snow.jpg')
plt.imshow(image, cmap='gray')
plt.gcf().set_size_inches(8, 8)
No description has been provided for this image
In [3]:
# Corrupt the image with Gaussian noise (provided)
image_noisy = noise.add_noise(image, 'gaussian')
plt.imshow(image_noisy, cmap='gray')
plt.gcf().set_size_inches(8, 8)
No description has been provided for this image

Note: from now on, please use the noisy image as the input for the filters.¶

1.1 Filter the noisy image with a 3x3 moving average filter. Show the filtering results.¶

In [4]:
# Design the filter h
### Insert your code ###
h = np.full((3, 3), 1/9)

# Convolve the corrupted image with h using scipy.signal.convolve2d function
### Insert your code ###
image_filtered = scipy.signal.convolve2d(image_noisy, h, mode='same')

# Print the filter (provided)
print('Filter h:')
print(h)

# Display the filtering result (provided)
plt.imshow(image_filtered, cmap='gray')
plt.gcf().set_size_inches(8, 8)

Filter h:
[[0.11111111 0.11111111 0.11111111]
 [0.11111111 0.11111111 0.11111111]
 [0.11111111 0.11111111 0.11111111]]
No description has been provided for this image

1.2 Filter the noisy image with a 11x11 moving average filter.¶

In [5]:
# Design the filter h
### Insert your code ###
h = np.full((11, 11), 1/121)

# Convolve the corrupted image with h using scipy.signal.convolve2d function
### Insert your code ###
image_filtered = scipy.signal.convolve2d(image_noisy, h, mode='same')

# Print the filter (provided)
print('Filter h:')
print(h)

# Display the filtering result (provided)
plt.imshow(image_filtered, cmap='gray')
plt.gcf().set_size_inches(8, 8)

Filter h:
[[0.00826446 0.00826446 0.00826446 0.00826446 0.00826446 0.00826446
  0.00826446 0.00826446 0.00826446 0.00826446 0.00826446]
 [0.00826446 0.00826446 0.00826446 0.00826446 0.00826446 0.00826446
  0.00826446 0.00826446 0.00826446 0.00826446 0.00826446]
 [0.00826446 0.00826446 0.00826446 0.00826446 0.00826446 0.00826446
  0.00826446 0.00826446 0.00826446 0.00826446 0.00826446]
 [0.00826446 0.00826446 0.00826446 0.00826446 0.00826446 0.00826446
  0.00826446 0.00826446 0.00826446 0.00826446 0.00826446]
 [0.00826446 0.00826446 0.00826446 0.00826446 0.00826446 0.00826446
  0.00826446 0.00826446 0.00826446 0.00826446 0.00826446]
 [0.00826446 0.00826446 0.00826446 0.00826446 0.00826446 0.00826446
  0.00826446 0.00826446 0.00826446 0.00826446 0.00826446]
 [0.00826446 0.00826446 0.00826446 0.00826446 0.00826446 0.00826446
  0.00826446 0.00826446 0.00826446 0.00826446 0.00826446]
 [0.00826446 0.00826446 0.00826446 0.00826446 0.00826446 0.00826446
  0.00826446 0.00826446 0.00826446 0.00826446 0.00826446]
 [0.00826446 0.00826446 0.00826446 0.00826446 0.00826446 0.00826446
  0.00826446 0.00826446 0.00826446 0.00826446 0.00826446]
 [0.00826446 0.00826446 0.00826446 0.00826446 0.00826446 0.00826446
  0.00826446 0.00826446 0.00826446 0.00826446 0.00826446]
 [0.00826446 0.00826446 0.00826446 0.00826446 0.00826446 0.00826446
  0.00826446 0.00826446 0.00826446 0.00826446 0.00826446]]
No description has been provided for this image

1.3 Comment on the filtering results. How do different kernel sizes influence the filtering results?¶

Insert your answer¶

A larger kernel is better at removing noise because it averages over more pixels, using more information to better estimate the true value of each pixel. When averaging, a larger dataset reduces the impact of outliers (noisy pixels, in this context) on the final result.

However, a larger kernel also significantly blurs the image. This happens because averaging over a wider area causes pixel values to converge toward the average of the surrounding pixels, reducing fine detail.

Averaging involves a trade-off between noise reduction and detail preservation. Larger kernels remove more noise but sacrifice more detail. A more effective alternative might be a median filter with a small kernel size, which removes noise efficiently (being less sensitive to outliers) while preserving detail.

Additionally, larger kernels are computationally slower, as they involve averaging over a greater number of pixels, increasing the computational workload.

2. Edge detection (56 points).¶

Perform edge detection using Sobel filtering, as well as Gaussian + Sobel filtering.

2.1 Implement 3x3 Sobel filters and convolve with the noisy image.¶

In [7]:
# Design the filters
### Insert your code ###
sobel_x = np.array([[1, 0, -1], [2, 0, -2], [1, 0, -1]])
sobel_y = np.array([[1, 2, 1], [0, 0, 0], [-1, -2, -1]])

# Image filtering
### Insert your code ###
g_x = scipy.signal.convolve2d(image_noisy, sobel_x, mode='same')
g_y = scipy.signal.convolve2d(image_noisy, sobel_y, mode='same')

# Calculate the gradient magnitude
### Insert your code ###
grad_mag = np.sqrt(g_x**2 + g_y**2)

# Print the filters (provided)
print('sobel_x:')
print(sobel_x)
print('sobel_y:')
print(sobel_y)

# Display the magnitude map (provided)
plt.imshow(grad_mag, cmap='gray')
plt.gcf().set_size_inches(8, 8)

sobel_x:
[[ 1  0 -1]
 [ 2  0 -2]
 [ 1  0 -1]]
sobel_y:
[[ 1  2  1]
 [ 0  0  0]
 [-1 -2 -1]]
No description has been provided for this image

2.2 Implement a function that generates a 2D Gaussian filter given the parameter $\sigma$.¶

In [8]:
# Design the Gaussian filter
def gaussian_filter_2d(sigma):
    # sigma: the parameter sigma in the Gaussian kernel (unit: pixel)
    #
    # return: a 2D array for the Gaussian kernel
    
    ### Insert your code ###
    # ignoring small values (outside k*sigma away from centre)
    k = 3
    
    # create grid of (i, j) coords
    x = np.arange(-k*sigma, k*sigma+1)
    y = np.arange(-k*sigma, k*sigma+1)
    i, j = np.meshgrid(x, y)

    # Calculate the Gaussian function
    h = 1/(2*math.pi*sigma**2) * np.exp(-(i**2 + j**2) / (2*sigma**2))
    
    return h

# Visualise the Gaussian filter when sigma = 5 pixel (provided)
sigma = 5
h = gaussian_filter_2d(sigma)
plt.imshow(h)
Out[8]:
<matplotlib.image.AxesImage at 0x24cd7c07c50>
No description has been provided for this image

2.3 Perform Gaussian smoothing ($\sigma$ = 5 pixels) and evaluate the computational time for Gaussian smoothing. After that, perform Sobel filtering and show the gradient magintude map.¶

In [9]:
# Construct the Gaussian filter
### Insert your code ###
h = gaussian_filter_2d(5)

# Perform Gaussian smoothing and count time
### Insert your code ###
start_time = time.time()
image_smooth = scipy.signal.convolve2d(image_noisy, h, mode='same')
end_time = time.time()
print(f'Computational time for Gaussian smoothing: {end_time - start_time} seconds')

# Image filtering
### Insert your code ###
g_x = scipy.signal.convolve2d(image_smooth, sobel_x, mode='same')
g_y = scipy.signal.convolve2d(image_smooth, sobel_y, mode='same')

# Calculate the gradient magnitude
### Insert your code ###
grad_mag = np.sqrt(g_x**2 + g_y**2)

# Display the gradient magnitude map (provided)
plt.imshow(grad_mag, cmap='gray', vmin=0, vmax=100)
plt.gcf().set_size_inches(8, 8)

Computational time for Gaussian smoothing: 2.756606101989746 seconds
No description has been provided for this image

2.4 Implement a function that generates a 1D Gaussian filter given the parameter $\sigma$. Generate 1D Gaussian filters along x-axis and y-axis respectively.¶

In [10]:
# Design the Gaussian filter
def gaussian_filter_1d(sigma):
    # sigma: the parameter sigma in the Gaussian kernel (unit: pixel)
    #
    # return: a 1D array for the Gaussian kernel
    
    ### Insert your code ###
    k = 3
    x = np.arange(-k*sigma, k*sigma+1)
    h = 1/(np.sqrt(2*math.pi)*sigma) * np.exp(-(x**2) / (2*sigma**2))
    return h

# sigma = 5 pixel (provided)
sigma = 5

# The Gaussian filter along x-axis. Its shape is (1, sz).
### Insert your code ###
h_x = gaussian_filter_1d(sigma).reshape(1, -1)

# The Gaussian filter along y-axis. Its shape is (sz, 1).
### Insert your code ###
h_y = gaussian_filter_1d(sigma).reshape(-1, 1)

# Visualise the filters (provided)
plt.subplot(1, 2, 1)
plt.imshow(h_x)
plt.subplot(1, 2, 2)
plt.imshow(h_y)
Out[10]:
<matplotlib.image.AxesImage at 0x24cdcfbc4d0>
No description has been provided for this image

2.5 Perform Gaussian smoothing ($\sigma$ = 5 pixels) using two separable filters and evaluate the computational time for separable Gaussian filtering. After that, perform Sobel filtering, show the gradient magnitude map and check whether it is the same as the previous one without separable filtering.¶

In [11]:
# Perform separable Gaussian smoothing and count time
### Insert your code ###
start_time = time.time()
image_smooth_x = scipy.signal.convolve2d(image_noisy, h_x, mode='same')
image_smooth = scipy.signal.convolve2d(image_smooth_x, h_y, mode='same')
end_time = time.time()
print(f'Computational time for separable Gaussian smoothing: {end_time - start_time} seconds')

# Image filtering
### Insert your code ###
g_x = scipy.signal.convolve2d(image_smooth, sobel_x, mode='same')
g_y = scipy.signal.convolve2d(image_smooth, sobel_y, mode='same')

# Calculate the gradient magnitude
### Insert your code ###
grad_mag2 = np.sqrt(g_x**2 + g_y**2)

# Display the gradient magnitude map (provided)
plt.imshow(grad_mag2, cmap='gray', vmin=0, vmax=100)
plt.gcf().set_size_inches(8, 8)

# Check the difference between the current gradient magnitude map
# and the previous one produced without separable filtering. You
# can report the mean difference between the two.
### Insert your code ###
mean_diff = np.mean(np.abs(grad_mag - grad_mag2))
print(f'Mean difference between the two gradient magnitude maps: {mean_diff}')

Computational time for separable Gaussian smoothing: 0.35946154594421387 seconds
Mean difference between the two gradient magnitude maps: 4.246997790888335e-13
No description has been provided for this image

2.6 Comment on the Gaussian + Sobel filtering results and the computational time.¶

Insert your answer¶

Results:

  • Mean difference between gradient magnitude maps produced by non-separable and separable Gaussian + Sobel: $4.25×10^{−14}$
  • Computation time for non-separable Gaussian + Sobel: $2.66s$
  • Computation time for separable Gaussian + Sobel: $0.363s$

(Note: These numbers are averages over 10 runs performed on my (fs422's) laptop, rounded to three significant figures.)

The reported mean difference between the two methods is negligibly small, likely caused by floating-point errors. This indicates that the two methods produce effectively identical results, as expected since the Gaussian filter is mathematically separable.

The separable Gaussian filter is over 7x faster than the non-separable Gaussian filter. This performance difference arises because the separable filter applies two 1D convolutions, while the non-separable filter performs a single 2D convolution. The computational complexity of a 2D convolution is $O(n2)$ per pixel, while two 1D convolutions have a combined complexity of $O(2n)$, making the separable approach significantly more efficient, with greater efficiency gains as the kernel size increases.

3. Challenge: Implement 2D image filters using Pytorch (24 points).¶

Pytorch is a machine learning framework that supports filtering and convolution.

The Conv2D operator takes an input array of dimension NxC1xXxY, applies the filter and outputs an array of dimension NxC2xXxY. Here, since we only have one image with one colour channel, we will set N=1, C1=1 and C2=1. You can read the documentation of Conv2D for more detail.

In [12]:
# Import libaries (provided)
import torch

3.1 Expand the dimension of the noisy image into 1x1xXxY and convert it to a Pytorch tensor.¶

In [13]:
# Expand the dimension of the numpy array
### Insert your code ###
image_noisy_expanded = np.reshape(image_noisy, (1, 1) + image_noisy.shape)

# Convert to a Pytorch tensor using torch.from_numpy
### Insert your code ###
image_noisy = torch.from_numpy(image_noisy_expanded)

3.2 Create a Pytorch Conv2D filter, set its kernel to be a 2D Gaussian filter and perform filtering.¶

In [14]:
# A 2D Gaussian filter when sigma = 5 pixel (provided)
sigma = 5
h = gaussian_filter_2d(sigma)

# Create the Conv2D filter
### Insert your code ###
h_expanded = np.reshape(h, (1, 1) + h.shape)
guassian_layer = torch.nn.Conv2d(1, 1, kernel_size=h.shape, padding=h.shape[0]//2, bias=False)
guassian_layer.weight.data = torch.from_numpy(h_expanded)

# Filtering
### Insert your code ###
image_filtered = guassian_layer(image_noisy)

# Display the filtering result (provided)
plt.imshow(image_filtered.squeeze().detach().numpy(), cmap='gray')
plt.gcf().set_size_inches(8, 8)
No description has been provided for this image

3.3 Implement Pytorch Conv2D filters to perform Sobel filtering on Gaussian smoothed images, show the gradient magnitude map.¶

In [15]:
# Create Conv2D filters
### Insert your code ###
sobel_x_layer = torch.nn.Conv2d(1, 1, 3, 1, bias=False)
sobel_y_layer = torch.nn.Conv2d(1, 1, 3, 1, bias=False)

sobel_x_layer.weight.data = torch.from_numpy(sobel_x).float().unsqueeze(0).unsqueeze(0)
sobel_y_layer.weight.data = torch.from_numpy(sobel_y).float().unsqueeze(0).unsqueeze(0)

# Perform filtering
### Insert your code ###
image_filtered_float = image_filtered.float()
g_x = sobel_x_layer(image_filtered_float)
g_y = sobel_y_layer(image_filtered_float)

# Calculate the gradient magnitude map
### Insert your code ###
grad_mag3 = torch.sqrt(g_x**2 + g_y**2).squeeze().detach().numpy()

# Visualise the gradient magnitude map (provided)
plt.imshow(grad_mag3, cmap='gray', vmin=0, vmax=100)
plt.gcf().set_size_inches(8, 8)
No description has been provided for this image
In [ ]: